Kelly_formula

#find x to maximized f
#f = (1+bx)^p(1-x)^q
#ifelse(g,bx,-x),P(g==T) = p, P(g==F) = q

n = 100
p = 0.7 # 
b = 0.8 #return from investing 1 unit, odds b to 1
kf = (b*p-(1-p))/b #Kelly formula, the fraction of total capital
if(kf < 0){
   kf = 0
   warning("warning: not to invest ")
}

x = 0.02
er = 1+x*(p*b+p-1)
#(1+b*x)*p+(1-x)*(1-p) 
#p+b*p*x+1-x-p+x*p = 1-x+p*x*b+p*x

vr2 = x^2*(1+b)^2*p*(1-p)
#(1+b*x - 1+x-p*x*b-p*x)^2 * p +(1-x - 1+x-p*x*b-p*x)^2*(1-p) 
#((1-p)*x*(b+1))^2*p + (p*x*(1+b))^2*(1-p)

sdr = x*(1+b)*sqrt(p*(1-p)) # linear to x

assumed you invest 20% of your capital. then after one period, your expected capital is 1.052 with standard deviation 0.1649727

#the function kelly try to optimized -- the accumulated return.
# as the number of bets increase, W/N -> p, and L/N -> 1-p. so it converge to this function
zf = function(r,b=1.6,p=0.6){(1+b*r)^p*(1-r)^(1-p)}


#in term of the tree method, for n bets
zf.tree = function(x,b=1.6,p=0.6,n=100){
   ss = c(1+b*x,1-x)
   ps = dbinom(0:n,n,p)
   A = c(1,rep(ss[1],n)) %>% cumprod()
   B = c(1,rep(ss[2],n)) %>% cumprod() %>% rev()
   return(list(ps = ps, values = A*B))
}


#in term of the tree method, for n bets
zft = function(x,b=1.6,p=0.6,n=100){
   ans = zf.tree(x,b,p,n)
   ps = ans$ps
   Vs = ans$values
   e.v = sum(ps*Vs)
   e.s2 = sum(ps*(Vs-e.v)^2)
   return(list(e.ret = e.v, e.sd = sqrt(e.s2)))
}

#in term of the tree method, for n bets
zft.log = function(x,b=1.6,p=0.6,n=100){
   ans = zf.tree(x,b,p,n)
   ps = ans$ps
   Vs = log(ans$values) #use log transform, because the upper end is too big. take it as utility function
   e.v = sum(ps*Vs)
   e.s2 = sum(ps*(Vs-e.v)^2)
   return(list(e.ret = e.v, e.sd = sqrt(e.s2)))
}

zft.mr = function(x,b=1.6,p=0.6,n=100,a=0.95){
   ans = zf.tree(x,b,p,n)
   psc = cumsum(ans$ps)
   e.v = ans$values[which.max(ans$ps)]
   a1 = (1-a)/2
   a2 = 1-a1
   ind1 = which.min(abs(psc -a1))
   ind2 = which.min(abs(psc -a2))
   e.sd = ans$value[ind2]-ans$value[ind1]
   return(list(e.ret = e.v, e.sd = e.sd, e.sd_low = ans$value[ind1], e.sd_hi = ans$value[ind2]))
}

zft.mm = function(x,b=1.6,p=0.6,n=100,w=0.5){
   ans = zf.tree(x,b,p,n)
   ind = which.max(ans$ps)
   e.v = ans$values[ind]
   c1 = 0
   ind1 = ind-1
   while((ind1>0) & (c1 < w/2)){
      c1 = c1+ans$ps[ind1]
      ind1 = ind1-1
   }
   c2 = 0
   ind2 = ind+1
   while((ind2<n+1) & (c2 < w/2)){
      c2 = c2+ans$ps[ind2]
      ind2 = ind2+1
   }
   if(ind1 == 0) ind1 = ind1+1
   if(ind2 >n) ind2 = ind2-1
   e.sd = ans$value[ind2]-ans$value[ind1]

   #print(c(ind,ans$ps[ind],ind1,ind2,ans$ps[ind1:ind2]))
   return(list(e.ret = e.v, e.sd = e.sd, e.sd_low = ans$value[ind1], e.sd_hi = ans$value[ind2]))
}

#x is the fraction to be used for the bet
gdx_data = function(p,b,x,m_paths=500,n_steps = 100){
   ps = c(p,1-p)
   ss = c(1+b*x,1-x)
   r.ind = sample(c(1,2),m_paths*n_steps,replace = T, prob = ps)
   r = ss[r.ind]
   gd = data.frame(path = gl(m_paths,n_steps),d.ret = r,o = rep(1:n_steps,m_paths))
   gd = gd %>% group_by(path) %>% mutate(c.ret = cumprod(d.ret))
   gdx = rbind(data.frame(path = gl(m_paths,1),
                          d.ret = rep(1,m_paths),
                          o = rep(0,m_paths),
                          c.ret = rep(1,m_paths)),gd)
   return(gdx)
}

data = gdx_data(p,b,kf,m_paths = 800)
data$x = kf
xs = c(0.02,0.2,0.5,0.7,0.8)
for(x in xs){
   gdx = gdx_data(p,b,x)
   gdx$x = x
   data = rbind(data,gdx)
}

gp = ggplot(data = data,aes(x = o, y = c.ret, group = path))  + 
      geom_line(alpha = 0.02) +scale_y_log10()+facet_wrap(~x)
print(gp)

the theoratical distribution of different fractions to invest

• as the invested amount increasing, the variance seems to be increasing.but it’s log scale, the lower side is near zeros. if the range to be considered is narrow, the up side would not include the extreme cases.
• but the maximum likely outcome, or the peak shifts to right but falls back to the left after a critical value.
• the k for maximum likely outcome can be found in range: \[np - q \leq k \leq np+p \]
• to track the peak can be done by k = np
• the maximum likely value is: \[(1+bx)^{np}(1-x)^{n-np} = [(1+bx)^p(1-x)^q]^n\], that is the Kelly function
• so, Kelly’s formula maximized the peak position.

x = 0.2
gf = function(x){zf.tree(x,b,p) %>% as.data.frame() %>% mutate(x = x, pc = cumsum(ps),cf95 = pc > 0.025 & pc <0.975 )}
g = gf(x)
ps = g$ps
ind = which.min(abs(cumsum(ps)-0.975))
p.level = ps[ind]
xs = c(0.02,0.2,0.5,0.7,0.8,kf)
g.data = map_dfr(xs,gf)
pg = ggplot(g.data %>% filter(values > 1e-15 & values < 1e6),aes(x = values, y = ps,color = factor(x)))+
   geom_point(alpha = 0.7,size = 0.5)+geom_line(aes(group = factor(x)),alpha = 0.5)+
   geom_vline(xintercept = 1,linetype="dashed")+
   geom_hline(yintercept = p.level,color="red",alpha = 0.4,linetype="dashed")+
   #facet_wrap(~factor(x),ncol = 1)+
   scale_x_log10()+theme(legend.position = "none")
print(pg)

a few points

• the normal average is dominated by the absurdly large number in the upper side, even the probability of happening is so slim.
• the log transform or the utility function lower the weights of those large numbers, and bring the expectation reflects more about the values in the common range.
• in the table below, the log-transformed mean is more close to the median values.
• there are simulation data and theoretical data.

data.end = data %>% filter(o == max(o))  %>% mutate(log.c.ret = log(c.ret))
data.end.stat = data.end %>% group_by(x) %>% 
      summarise(e.ret = mean(c.ret),s.ret = sd(c.ret),
                m.ret = quantile(c.ret,0.5),
                q05.ret = quantile(c.ret,0.05),
                q95.ret = quantile(c.ret,0.95),
                e.elog.ret = mean(log.c.ret)%>%exp(),s.elog.ret = sd(log.c.ret)%>%exp(),
                .groups = "keep")
knitr::kable(data.end.stat,digits = 3)

x	e.ret	s.ret	m.ret	q05.ret	q95.ret	e.elog.ret	s.elog.ret
0.020	1.679	0.278	1.657	1.242	2.137	1.656	1.180
0.200	129.715	288.974	40.249	2.940	786.508	36.624	5.153
0.325	1491.265	6923.664	80.362	0.545	6346.273	77.107	15.409
0.500	21932.848	122078.999	15.778	0.011	59611.065	21.489	121.295
0.700	2614.954	32623.556	0.007	0.000	698.848	0.003	2268.210
0.800	666.551	8393.498	0.000	0.000	2.929	0.000	20291.095

gp = data.end %>%
      ggplot(aes(x = x, y = c.ret))  + 
      stat_summary(alpha = 0.8,geom = "bar",fun = "mean",width = 0.03)+
      stat_summary(geom = "errorbar",fun.data = "mean_se",width = 0.01)+
      geom_hline(yintercept = 1,linetype = "dashed")+
      scale_y_log10()
print(gp)

#p = 0.7 # 
#b = 0.7 #return from investing 1 unit, odds b to 1
#er = b*p-1*(1-p) #expected return
#kf = (b*p-(1-p))/b 

xs = c(0.02,0.2,0.5,0.7,0.8)
zx = seq(0,0.95,0.01)
zd = data.frame(fraction = zx, expected.return = zf(zx,b,p),sdr = zx*(1+b)*sqrt(p*(1-p))/10)
zd2 = data.frame(fraction = xs, expected.return = zf(xs,b,p),sdr = xs*(1+b)*sqrt(p*(1-p))/10)
zd$color = "grey"
zd2$color = "blue"
zd = rbind(zd,zd2)

pg = ggplot(data = zd)+geom_point(aes(x = fraction, y = expected.return,color = color))+
         geom_line(aes(x = fraction, y = expected.return+sdr),color = "grey")+
         geom_line(aes(x = fraction, y = expected.return-sdr),color = "grey")+
         geom_point(aes(x = kf,y = zf(kf,b,p)),color = "red")+
         scale_colour_identity()
print(pg)

#p = 0.7 # 
#b = 0.7 #return from investing 1 unit, odds b to 1
#er = b*p-1*(1-p) #expected return
#kf = (b*p-(1-p))/b 

xs = c(0.02,0.2,0.5,0.7,0.8)
zx = seq(0,0.8,0.01)
f.zft = function(x){zft.log(x,b,p)  }

zd =  map_df(zx,f.zft) %>% rename(expected.return = e.ret, sdr = e.sd) %>% 
         mutate(fraction = zx,sdr = sdr,color = "grey")
zd2 =  map_df(xs,f.zft) %>% rename(expected.return = e.ret, sdr = e.sd) %>% 
         mutate(fraction = xs,sdr = sdr,color = "blue")

zd = rbind(zd,zd2)

pg = ggplot(data = zd)+geom_point(aes(x = fraction, y = expected.return,color = color))+
         geom_line(aes(x = fraction, y = expected.return+sdr),color = "grey")+
         geom_line(aes(x = fraction, y = expected.return-sdr),color = "grey")+
         geom_point(aes(x = kf,y = f.zft(kf)$e.ret),color = "red")+
         scale_colour_identity()+ylab("expected.log.return")
print(pg)

pg = ggplot(data = zd)+geom_point(aes(x = fraction, y = expected.return %>% exp(),color = color))+
         geom_line(aes(x = fraction, y = (expected.return+sdr/10)%>% exp()),color = "grey")+
         geom_line(aes(x = fraction, y = (expected.return-sdr/10)%>% exp()),color = "grey")+
         geom_point(aes(x = kf,y = f.zft(kf)$e.ret %>% exp()),color = "red")+
         scale_colour_identity()+ylab("expected.log.return")
print(pg)

#p = 0.7 # 
#b = 0.7 #return from investing 1 unit, odds b to 1
#kf = (b*p-(1-p))/b 
#if(kf < 0) kf = 0
xs = c(0.02,0.2,0.5,0.7,0.8)
zx = seq(0,0.8,0.01)
f.zft = function(x){zft.mr(x,b,p,a = 0.3)} #the central 30% confident zone

zd =  map_df(zx,f.zft) %>% rename(expected.return = e.ret, sdr = e.sd) %>% 
         mutate(fraction = zx,sdr = log1p(sdr),color = "grey")
zd2 =  map_df(xs,f.zft) %>% rename(expected.return = e.ret, sdr = e.sd) %>% 
         mutate(fraction = xs,sdr = log1p(sdr),color = "blue")

zd = rbind(zd,zd2)

pg = ggplot(data = zd)+geom_point(aes(x = fraction, y = expected.return,color = color))+
         geom_line(aes(x = fraction, y = e.sd_low),color = "grey")+
         geom_line(aes(x = fraction, y = e.sd_hi),color = "grey")+
         geom_point(aes(x = kf,y = f.zft(kf)$e.ret),color = "red")+
         scale_colour_identity()+ylab("most.likely.return")


print(pg+scale_y_log10())

print(pg)

the relationship between variance and return

• in the normal mean calculation, the mean and standard devision is positive correlated
• in the log mean calculation, there is optimal point for maximizing the return while minimizing the variance.
• in the confident range method, if the range is too narrow the extreme case would not take into account. the difference between the upper value and lower value would much smaller.

zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft(x,b,p)}
d = map_df(zx,f.zft) %>% mutate(fraction = zx)
pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
      scale_y_log10()+scale_x_log10()
print(pg)

zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft.log(x,b,p)}
d = map_df(zx,f.zft) %>% mutate(fraction = zx)
pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
      geom_point(aes(x = f.zft(kf)$e.sd, y = f.zft(kf)$e.ret),color = "red")

print(pg)

zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft.mr(x,b,p,a = 0.95)}
d = map_df(zx,f.zft) %>% mutate(fraction = zx)


pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
      geom_point(aes(x = f.zft(kf)$e.sd, y = f.zft(kf)$e.ret),color = "red")

print(pg)

zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft.mr(x,b,p,a = 0.3)}
d1 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d1$a = 0.3

f.zft = function(x){zft.mr(x,b,p,a = 0.6)}
d2 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d2$a = 0.6

f.zft = function(x){zft.mr(x,b,p,a = 0.1)}
d3 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d3$a = 0.1

d = rbind(d1,d2,d3)

pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
      geom_path(aes(x = e.sd,y = e.ret,linetype = factor(a)),color = "grey")

print(pg)

zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft.mm(x,b,p,w = 0.9)}
d = map_df(zx,f.zft) %>% mutate(fraction = zx)


pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
      geom_point(aes(x = f.zft(kf)$e.sd, y = f.zft(kf)$e.ret),color = "red")

print(pg)

zx = seq(0.01,0.8,0.01)
f.zft = function(x){zft.mm(x,b,p,w = 0.3)}
d1 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d1$w = 0.3

f.zft = function(x){zft.mm(x,b,p,w = 0.6)}
d2 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d2$w = 0.6

f.zft = function(x){zft.mm(x,b,p,w = 0.1)}
d3 = map_df(zx,f.zft) %>% mutate(fraction = zx)
d3$w = 0.1

d = rbind(d1,d2,d3)

pg = ggplot(d)+geom_point(aes(x = e.sd,y = e.ret,color = fraction))+
      geom_path(aes(x = e.sd,y = e.ret,linetype = factor(w)),color = "grey")

print(pg)

the fair bet

• the direct draw up for the kelly function is from 1 to 1+bx and 1-x with probability p and 1-p.
• the bet for kelly function is bet x on the b to 1 against bet. that is x to (1+b)x and 0 with probability p and 1-p.
• fair bet condition: \[(1+b)*p = 1\]
• \(p\leq1\) and \(b \gt 0\) so, \[k_f = p - \frac{(1-p)}{b} \leq1\]

DiagrammeR::grViz("digraph {
  node [shape = oval]
  a [label = '1']
  b [label = '1+bx']
  c [label = '1-x']
  a->b [label='p']
  a->c [label='1-p']
  
  d [label = 'x']
  e [label = '(1+b)x']
  f [label = '0']
  d->e [label='p']
  d->f [label='1-p']
  {rank=same;a;d;}
  }",
  height = 200)

N = 1e3
#pb.data = data.frame(p = runif(N),b = rnorm(N,mean = 1,sd = 0.3)) %>% filter(b >0)
pb.data = data.frame(p = runif(N),b = runif(N,min = 0,max = 2))%>% filter(b >0)

pb.data = pb.data %>% mutate(kf = (b*p-(1-p))/b,jd = kf > 0,kf_ = ifelse(kf>0,kf,0))
pg = ggplot(pb.data)+geom_point(aes(x = p, y = b,shape = jd,color = kf_))+
      geom_line(data = data.frame(x = seq(0.33,1,0.01)),aes(x = x,y = 1/x-1),color = "red")+
      guides(shape=guide_legend(title="kf > 0"))+
      scale_shape_manual(values=c(1,3))
print(pg)